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Abstract 

The long time dynamics of polymeric materials has been extensively studied in the past through various 
experimental techniques and computer simulations. While computer simulations typically treat generic, 
simplified models, experiments deal with specific chemistries. In this letter we present a hierarchical ap- 
proach that combines atomistic and coarse-grained simulations to quantitatively study polymer dynamics. 
As an example we predict diffusion coefficients of atactic polystyrene melts of molecular weights relevant to 
polymer processing (up to 50kDa) without any adjustable parameter and compare the results to experiment. 

PACS numbers: 36.20-y, 61.25.H-, 83.10Rs 
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Introduction.— The dynamics of polymers is a prototypical soft matter prob 



em and has at- 



tracted considerable experimental and theoretical attention for many years yj, |2J, [3|, |4j]. For the 
two limiting cases of short, but already fully flexible polymer chains in a melt, and very long, 
fully entangled chains well established (scaling) theories, namely the Rouse and reptation theory, 
exist and have been tested in detail by experiments and computer simulations of highly idealized 



models il, Q, El E]. However most of these approaches provide either a semi quantitative link 
between theory and experiment in the sense that scaling relations are verified, or quantitative in- 
formation for short chains with simple chemical structure, like PE. What is still missing is a direct 
quantitative link between chemical structure and measurable quantities like diffusion constants 
or viscosities, especially for high molecular weight entangled systems. In addition the crossover 
regime from Rouse like to reptation like behavior, which for example for polystyrene (PS) covers 
a range of molecular weights from about M RilkDa to M ^lOOkDa, is still under discussion. 
In this regime it is also not clear, to what extent universal theories describe the properties and to 
what extent chemical details surface strongly, which makes the applicability of a general analytic 
theory virtually impossible. This is a situation, where well tailored computer simulation can be of 
significant help. 

On the microscopic level detailed atomistic molecular dynamics (MD) simulations allow quan- 
titative predictions of the dynamics of simple polymers with rather low molecular weight [(]]. 
However, due to the broad spectrum of characteristic lengths and times involved in polymeric sys- 
tems, it is not feasible to apply them to systems of more complex chemical structure or of high 
molecular weight. 

On the mesoscopic level coarse-grained (CG) molecular dynamics and Monte Carlo simula- 
tions have proven to be very efficient means to study the dynamics of long, entangled simple 
model polymer systems [|5|,LZ|]. In a similar way, structure-based CG models have been employed 
to study the dynamics of specific polymer melts. There groups of chemically connected atoms are 
lumbed into 'superatoms'. Effective, coarse-grained interaction potentials are obtained by averag- 
ing over microscopic details of the underlying atomistic models. The direct link to the chemistry 
is maintained through the choice of appropriate superatoms and the resulting set of structurally de- 
fined CG bonded and non-bonded temperature dependent effective potentials (more precisely free 
energies) [5,8,^]. By doing that structural properties of polymeric systems are described quite 
well. However CG simulations cannot be used for a direct quantitative study of dynamics because 
the intrinsic time scale of the CG model is not the same as that of the underlying chemical system. 
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One reason is the use of (softer) effective CG potentials, which results in a significantly reduced 
effective friction between the beads compared to the original monomers. To overcome this lim- 
itation the ratio between the CG time and chemically realistic time scale has to be derived from 
either experiments or atomistic simulations [8|,I9D. However a direct comparison of the resulting 
dynamics to experiment and the extension to long, entangled chains are still missing. 

In this letter we extend this approach and present a hierarchical methodology that combines 
dynamic simulations, on different length and time scales. This approach allows us to predict 
through CG simulations, not only qualitatively, but also quantitatively, the dynamics of polymer 
systems of molecular weight relevant to polymer processing without any adjustable parameter. As 
an example the whole methodology is applied in atactic polystyrene at T=463K, which has been 
extensively studied in the past through experiments (see for example [2L11O0 and references within 
and simulations [si, 9]). The chosen temperature is a typical process for PS. 

Simulation Methodology. — In order to derive a time scaling factor S between different models 
we match the mean square displacements of the monomers in amplitude and slope. The latter is 
important, since the motion characteristics of the different models coincides only above a charac- 
teristic length scale. A three level ansatz is needed in order to also capture the effect of the chain 
length dependent melt density onto the scaling factor. 

Detailed all atom simulations, using a model in which hydrogens and carbons are treated ex- 
plicitly II 1 lh . have been performed for rather short PS chains (M=lkDa, T=463K, only 10 repeat 
units) at the experimental density of 0.925gr/cm 3 . Only for such small chains reliable data for 
the mean square displacement of the monomers can be obtained with all-atom MD simulations. 
To somewhat overcome thisproblem united atom (UA) models are frequently employed; here we 
are using the TraPPE one lll2n . There hydrogens are lumped together with carbons, defining new 
CH Z types of united atoms. For UA simulations it is usually assumed that the time scale does not 
in a detectable way deviate from all atom simulations, because hydrogens are very small and the 
"coarse-graining" of the UA models is of the order of hydrogen-carbon bond which is only about 
1 A. However this assumption, which in many cases (for example UA models for PE, PB, PI) 
works reasonably well, fails to work for the UA-TraPPE PS model; i.e. it predicts for PS a much 
faster dynamics. This is most probably, due to the too low dihedral barriers or due the missing 
electrostatic interactions but does not affect the overall conformations [O]. Various systems, with 
molecular weight from lkDa to lOkDa (T=463K), have been studied by this method. Beyond that 
CG models are used, which however do not reproduce the motion patterns down to length scales 
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FIG. 1: Coarse-graining mapping scheme of PS: one monomer is mapped to two different CG beads 
{a a =4.1 A, mA=27amu and as =5. 2A, mB=77amu) |9|]. 



of the order of one or two monomers. 

Thus we adopt the following strategy. First UA and CG simulations are performed, for molecu- 
lar weights between 1 and lOkDa in order to obtain the time scaling factor Sjja-cg(M). Then the 
UA results are calibrated by all atom simulations, however for much shorter times only, resulting 
in an additional scaling factor S AA _ UA (M). With S(M) = S AA _ UA (M) ■ S UA _ C G(M)[sec/r] we 
calibrate the time scale to determine long chain polymer diffusion constants and compare these to 
experiment. 

The molecular dynamics package GROMACS 1130 was used to perform all the atomistic (AA 
and UA) MD simulations. Initial well-equilibrated atomistic polymer melts are obtained by back- 
mapping CG melts Q9[|. All simulations have been performed at constant temperature and volume 
(NVT ensemble) at experimental densities II 1-411 . The overall simulation time of the production 
runs ranged from 50 ns to 300 ns depending on the molecular weight. 

CG MD simulations have been performed using a CG model for PS in which one PS monomer 
is mapped onto two effective coarse grained beads, i.e. a 2:1 model 19fl (see Figured]). CG bead 
"A" corresponds to the CH 2 of a PS monomer plus the half mass of each of the two neighboring 
CH groups along the chain backbone, whereas CG bead "B" is the phenyl ring. This model 
does not loose too many structural details in comparison to all atom systems, while still being 
very efficient compared to atomistic simulations. Various monodisperse atactic PS melts, with 
molecular weight from lkDa up to 50kDa have been studied with CG MD simulations. Note 
that the characteristic molecular weight M e for the formation of entanglements for atactic PS is 

n 

about 15kDa (at T=463K) IIP)] . All CG systems have been originally generated by a combination 
of Monte Carlo and MD simulations following the procedure by Auhl et al. la, CG MD 



simulations are performed in dimensionless LJ units using m A to scale all masses, a AV = (a A + 



<7b)/2 to scale all lengths and e = kT to scale all energies. In order to control the temperature in 
the system we use a Langevin thermostat with friction coefficient T = l.Or -1 . All coarse-grained 



MD simulations are performed using the ESPResSO package [|16|l with a time step A?=0.01 r, 



with r = ^mAC^y/e. MD simulations ran for times between 1 x 10 4 r and 3 x 10 6 r depending 
upon the system size. 

Time Mapping. — The problem of time scales can be discussed in terms of local (monomeric) 
friction coefficients. For both the Rouse as well as the reptation model the local motion is governed 
by a scalar friction coefficient £, so that the melt viscosity 77 oc ( and the chain diffusion constant 
D(N) oc For the modeling of the polymer melt this bead friction depends on the specific 
representation of the polymer. The softer CG potentials result in a significantly reduced effective 
friction coefficient, £ CG , between beads compared to the friction coefficient in the united atom 
description, ( UA , and the all atom description, ( AA , respectively. As a consequence the time in the 
dynamic CG simulations does not correspond to the real time of the polymeric system and has to 
be properly scaled. In most cases the differences between all atom and united atom simulations are 
much smaller then the typical error bars of a simulation. For the TraPPE-UA PS model however 
this is not the case. Though in general a disadvantage, for PS this makes it possible to obtain a 
proper scaling by employing the AA-UA-CG hierarchy. Thus the time scaling parameter, which 
is the ratio of the effective bead frictions in the atomistic and the CG description, is given by 

S(M) = C AA /C CG = S A a-ua(M) ■ S UA _ CG (M) (1) 

Because the local energy landscape is quite complex and fluctuating it is not possible to give a 
well founded analytical prediction of S. Therefore, S should be obtained using data taken either 
from atomistic simulations, as done here, or directly from experiment. 

To determine S we match mean square displacements (MSD) of chain beads over a consider- 
able time, where both amplitude and shape coincide. First we examine the time mapping of the 
CG data based on the UA MD simulations. Figure [2] shows the mean square displacements aver- 
aged over all beads i of the CG model and the correspondingly analyzed UA model (circles), gi(t), 
gi (t) = {((rj(i) -ri(O)) 2 )}^) fromUA MD and CG MD simulations for a specific PS melt (lkDa, 
T=463K). The scaling factor, Sua-cg, m order to match the two curves on top of each other in the 
long time regime, gives SuA-CG(M=lkD&)= 3. lps/r. With this scaling factor both curves coincide 
above a distance of about 10 A and a corresponding time of about 200 ps. Below that distance and 
time the coarse graining results in a different shape of the curve, illustrating that a mere crossing 
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FIG. 2: Time scaling of the CG simulations using UA data and of the UA simulations using atomistic data 
for a PS melt (M=lkDa, T=463K), based on the motion of the polymer beads. 



of the curves is not sufficient to determine S. Note also that Sua-cg is the same if the mapping 
is based on a quantity describing the orientational dynamics (for example the end-to-end vector 
autocorrelation function) of the PS chains II 1711 . 

In general Sua-cg depends on molecular weight and density, which depends on M (chain end 
free volume effect), as shown in Figure [3] (circles) for the systems studied by both UA MD and 
CG MD simulations. As we can see Sua-cg varies in the short length regime (up to about 50 
monomers), ranging from 3.1ps/r to about 6.0ps/r, and then it remains constant. This is in phase 
with the observed change in density, which varies from 0.925g/cm 3 for lkDa to about 0.97g/cm 3 
for the lOkDa and higher molecular weight melts (see Figure [2]) [14]. At high molecular weights 
(above lOkDa) the change in the polymer dynamics is entirely due to the increase of the molecular 
weight. On the other hand, in the short length regime the density effect is very important. The 
latter one is not being described accurately in the CG model, resulting into a dependence of S on 
the density (and on the molecular length). 

The important result of Figure [3j is that a single value for the time scaling parameter S is 
appropriate to describe the dynamics of long polymer chains. However the UA model itself in- 
cludes some minor coarse graining and in the case of the TraPPE-UA PS model predicts faster 
dynamics than the all-atom simulations. Thus we have to follow a same procedure as above, how- 
ever now for the two models exhibiting atomistic detail. A typical example is shown in Figure 
[2] (squares). Though qualitatively similar, there is a remarkable quantitative difference. As Fig- 
ure [2] displays the two sets of data perfectly match from a distance above about 3A, this is for 
significantly smaller length scales than before, of the order of the size of a benzene ring (Rq is 
shown with arrow). For the present case of Figure |2]we arrive at Saa-u A(M=lkDa.)=35, resulting 
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FIG. 3: Time mapping of the CG simulations of the PS melts using UA and AA data, and density as a 
function of M (T=463K). 



in SAA-CG(M=YkDa)=l lOps/r. This together with computational efficiency of the TraPPE UA 
model allows us to determine the time scaling factor with its full chain length dependence (up to 
M=10kDa), as shown in Figure [3] (squares). Note the similar qualitative but the large quantitative 
difference between Sjja-cg(M) and Saa-cg(M). Alternatively one can follow the observation, 
that the variation of S follows approximately the changes in density rather than the molecular 
weight itself, even though this density change is due to the chain length variation. This is not sur- 
prising since, for polymers both density and monomeric friction coefficient vary with molecular 
weight because of the higher chain end free volume (free volume theory) and two melts with the 
same density would be expected to have the same friction coefficient [2] . Therefore by performing 
the time mapping for the short chain system but at the density of the longer chains one also can 
obtain a reliable estimate of Saa-ua- If we follow this procedure the combined time mapping 
Saa-cg(M) varies between ~ 1 lOps/r (for the lkDa system) and ~ 700ps/r for the high (lOkDa 
and above) molecular weight (polymeric) regime. Here we follow the latter because of the ex- 
tensive computational effort needed to obtain reliable data for MSDs of long all-atom PS chains. 
Saa-cg now can be applied to determine absolute diffusion constants in long chain PS melts, 
based on CG simulations. 

Self-Diffusion Coefficient of Polystyrene Melts. — Following the hierarchical methodology de- 
scribed above, the dynamical and rheological properties of long entangled polymeric systems can 
be obtained from the CG dynamic simulations in real units. Here we focus on the self-diffusion 
coefficient, D, which can be calculated directly from the MSD of the chain center-of-mass by 
rescaling r to sec by the above described procedure. Data for the self-diffusion coefficient of the 
CG PS melts, scaled with Saa-cg(M), as a function of the molecular weight (squares) are shown 
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FIG. 4: Self-diffusion coefficient of PS melts as a function of the molecular weight (T=463K). 



in Figure H] and are compared to experimental data (circles) Jlp|. Both simulation and experimen- 
tal data are not corrected for the chain end free volume. The range of molecular weights (up to 
50kDa) spans the regime from unentangled to entangled PS melts. 

The results show a remarkable qualitative and quantitative agreement between the experimental 
and the simulated diffusion coefficients (see figure H]). This is of particular importance if we con- 
sider that results from the CG dynamic simulations are compared to experimental data, by using 
only detailed atomistic simulations for a few reference short-chain systems, without any adjustable 
parameter. The error bars in the AA-scaled CG data are mainly due to the high error bars in the 
underlying all-atom data. The larger deviation between the simulation and the experimental data 
in the short length regime is not surprising if we consider the effect of the (small) polydispersity of 
the experimental data (/ ~ 1.04): in these short chains the presence of even only a small amount 
of PS oligomers acts like an effective dilution and can increase the free volume and the diffusion 
of the systems. 

Note also that the entire regime of both simulation and experimental data can be fitted using a 
power-law dependence (D ~ M~ b ) with a single exponent of about b ~ 2.1 ± 0.2. For short chain 
polymer melts the Rouse model predicts b=l while long entangled polymer melts the pure reptation 
theory predicts 6=2 yj]. However the Rouse model neglects the molecular weight dependence of 
the density and and thus of the friction coefficient. In order to eliminate this effect a correction 
can be made for D as it is often done in the analysis of experiments; this will be a part of a future 



work jn\\ . 



It is of interest to estimate the time scales involved in the dynamics of the PS systems studied 
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here: the diffusion coefficient of the higher molecular weight PS melts (50kDa) is of the order of 

about 10~ n cm 2 /sec. This results into a relaxation time of the whole chain - decorrelation of the 

ri 

end-to-end vector -, (according to reptation theory =< R 2 > /3n 2 D [1]) of about 6.0ms, 
many orders of magnitude longer that what can be modeled with atomistic molecular dynamic 
simulations in such systems. 

In conclusion we present a hierarchical simulation approach that combines dynamic simula- 
tions on different length and time scales. A mapping over a small range of molecular lengths, 
using atomistic and united-atom simulation data, shows that the time mapping parameter S varies 
as a function of chain length induced by different melt densities. The asymptotic plateau value 
of S can be used for scaling the CG dynamic results of longer polymeric chains, where it is not 
possible to have reliable atomistic data at all. 

Furthermore the multi-scale nature of the proposed scheme, combined with a proper back- 
mapping procedure BalSUla], allows to study time scales ranging from afewfs up to ms by MD 
simulations and to compare directly to experimental data without any adjustable parameter. The 
accuracy of the CG dynamical predictions, compared to experimental data, of course depends 
on the accuracy of the underlying microscopic (atomistic) model with which the CG "raw" data 
are scaled. This opens up the way for developing simulation methodologies that can be used for 
quantitative studies of the dynamics and the rheology of complex systems within the M regime 
relevant to the polymer processing. The proposed approach can be directly extended to other poly- 
mer chemistries by properly choosing a structure-based CG model, which reproduces structural 
properties for length scales of about 0.5- lnm and larger. Furthermore it can be used for the study 
such different systems as non-equilibrium polymer melts, polymers at temperatures near to T g , or 
polymer/solid interfacial systems. 

One the other hand, and despite the success of the proposed methodology, still a number of 
issues remain. For example, for multi-component systems or polymers in solution, the nature of 
the friction might not allow a renormalization in time with a single quantity. For such systems 
molecular dynamics simulations could be used for the calculation of the entire friction matrix, 
which can be then involved in the equations of motion on an even more coarse level, as for example 
in the level of the primitive paths IU9I1 or through the GENERIC formalism ll20n . 
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